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Abstract 

We study heat conduction in a billiard channel formed by two sinusoidal walls and the diffusion of 
particles in the corresponding channel of infinite length; the latter system has an infinite horizon, i.e, 
a particle can travel an arbitrary distance without colliding with the rippled walls. For small ripple 
amplitudes, the dynamics of the heat carriers is regular and analytical results for the temperature 
profile and heat flux are obtained using an effective potential. The study also proposes a formula 
for the temperature profile that is valid for any ripple amplitude. When the dynamics is regular, 
ballistic conductance and ballistic diffusion are present. The Poincare plots of the associated 
dynamical system (the infinitely long channel) exhibit the generic transition to chaos as ripple 
amplitude is increased. When no Kolmogorov-Arnold-Moser (KAM) curves are present to forbid 
the connection of all chaotic regions, the mean square displacement grows asymptotically with time 
t as t log (t) . 
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I. INTRODUCTION 



Much attention has been paid in the past decade to the study of heat conduction in ID 
and 2D systems via their underlying dynamics. The main question to be answered is how 
the Fourier law of heat conduction is related to the microscopic dynamics of the system. 
On the one hand, a convergent conductivity has been shown in the ding-a-ling model [l| 
(where oscillators exchange energy via intermediate hard spheres), which is chaotic, and in 
the Lorentz channel (a channel with circular scatterers placed periodically) which is fully 
hyperbolic j^j. On the other hand, for the Fermi-Pasta-Ulam chain, where oscillators are 
coupled to non-linear terms, heat conductivity is abnormal even above the chaotic threshold 

; while the serpent billiard, a channel with parallel semi-circular walls, exhibits marginally 
normal diffusion even under conditions of global chaos 4|. Since normal diffusion is at the 
root of normal heat transport and abnormal diffusion leads to abnormal heat transport jsj, 
the latter system cannot obey Fourier's law. Thus, the positivity of the Lyapunov exponent is 
neither a sufficient nor necessary condition for inducing normal transport properties, since 
there are billiard gas models (polygonal channels) with linear dynamical instability and 



yet they exhibit normal transport properties 6Hl0j]. In addition, it seems that for chaotic 



systems a strong degree of chaos, which translates into a global chaotic dynamics, is required 
to obtain normal transport properties. 

Although the precise conditions behind the onset of Fourier's law remain unknown despite 
decades of intensive studies, many interesting properties of the heat mechanism have been 
discovered, including the possibility of controlling heat flux: for example, there are ID 



chains that act as thermal rectifiers 



lJl3 



s where particle interactions (or 



13], billiard mode 

an external magnetic field) produce thermal rectification [l4,[l5| and graded systems in which 

3, IT 



rectification does not decay with system size \v\ . Moreover, a temporally alternatin. 
bath temperature can be used to generate a steady heat flow against a thermal bias 
while an important increasing of the thermoelectric efficiency is reported in billiard-like 
systems [19]. Another interesting, related issue is the possibility of controlling the stochastic 
transport of particles (in the absence of a thermal gradient) which can be achieved in systems 
possessing spatial or dynamical symmetry breaking with the aid of external unbiased input 



signals [20] 



Much attention has been paid to billiard systems when they exhibit strong chaos (global 



chaos), but few studies have focused on the relation between the degree of chaos in the sys- 



tem and their transport properties 
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22J, despite the importance of this issue, as in many 



systems chaotic and regular dynamics coexist. In addition, the system introduced in Ref. 



2l| is a peculiar chaotic system whose phase space does not have the typical KAM structure 



of generic Hamiltonian systems. While diffusion in systems (quasi-lD cosine billiard) with 
KAM phase space structure is studied in Ref. 22J, but no connection with the thermal 
transport properties is probed. Therefore, the main purpose of our study is to analyze the 
effect of the degree of chaos on the thermal transport properties of systems with KAM phase 
space structure. 

In this paper, we consider a two-dimensional billiard model formed by two sinusoidal walls. 
The two ends of the channel are connected to Gaussian-type thermal baths (see Fig. [TJ and 
the corresponding infinite length channel has an infinite horizon. The average width of the 
channel chosen is much smaller than its length in order to keep kinetic excitation in the 
transverse direction frozen. The study of thermal properties in this system is interesting 
because it exhibits Poincare sections with a KAM structure typical of generic Hamiltonian 
systems [23| . In addition, at small ripple amplitudes, we obtain an estimation of the temper- 
ature profile and heat flux by means of an effective potential. The result of our estimation 
is corroborated by numerical calculations. 

The paper is organized as follows. In Sec II we introduce the model and discuss the dy- 
namic properties of the system. Sec. Ill focuses on the study of the thermal transport 
properties of the channel when ripple amplitude is small. For this case, we use the effective 
potential mentioned in the previous section to obtain analytical results for the stationary 
heat flux and temperature profile. Section IV examines how the degree of chaos affects the 
thermal transport properties in the channel, focusing on the point at which the transition 
from regular dynamics to chaotic dynamics occurs. When the ripple amplitude is sufficiently 
large, the results are compared to those obtained for the Lorentz channel and for the system 
introduced in Ref. 21] . Finally, concluding remarks are presented. 
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II. SYSTEM MODEL AND DYNAMICS 

We considered an infinite channel formed by two sinusoidal walls; the profiles of the upper 
and lower walls, respectively, are determined by 

Hi = b + a sin(27rx) y 2 = — b + a sin(27r(x + r)), (1) 

where a is the amplitude of the ripples, r denotes the phase difference between the upper 
and lower walls, and b is half of the average width of the channel. All quantities are rescaled 
to the length I of one period. In this study, we consider the following fixed values: b = 0.1 
and r = 1/2. When we take other values for r, the main conclusions of this paper remain 
essentially the same. 

The dynamics of a particle that collides specularly with the rippled walls can be described 
qualitatively by the Poincare section (x n ,p n = cos(6' n )), where x n is the x-coordinate of the 
particle and 9 n the angle that the trajectory of the particle forms with the positive x-axis just 
after the n-th collision with any wall. To obtain all possible orbits in the Poincare section 
(PS), several initial conditions (x , yo, ) must be taken into account. Here, (x ,y ) is a 
point in the configuration space that corresponds to the initial x- and y-coordinates of the 
particle, while 6q is the angle that the initial trajectory of the particle forms with the positive 
x-axis at the departure point (x ,y ). Due to the periodicity of the channel, the structure 
of the PS is periodic; therefore, we may choose the x-domain in the interval [0,2], which 
corresponds to two ripple periods, so as to obtain a complete panorama of the dynamics 
of the system. Poincare plots for four different values of a are shown in Fig. [2j When 
the ripple amplitude is small, the dynamics of the system is regular and Poincare sections 
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resemble the phase space of a one- dimensional pendulum (see Fig. Et^a)). The elliptical 
orbits correspond to particles trapped in the channel that move adiabatically backwards 
and forwards around a stable fixed point. The trajectories outside the librational region 
represent a particle traveling to the left (p n < 0) or right (p n > 0) of the channel (rotational 
motion). When a is increased, the dynamics is still regular up to a < 0.015, and the region 
of librational motion occupies a larger size. For the interval 0.015 < a < 0.025 the separatrix 
becomes chaotic with some sizable width (see Fig. MJd)) that increases as ripple amplitude 
becomes larger. There are two KAM curves (KAM barriers) that forbid the connection of 
all chaotic regions. They are found in the limits between the red region (light gray) and 
the blue region (dark gray) shown in Fig|2](b). Therefore, motion is unidirectional for initial 
conditions lying outside both the librational region and the chaotic separatrix (we refer to 
this situation as unidirectional mixed chaos). There is a critical value of ripple amplitude 
(a m 0.025) at which all chaotic regions are connected due to the destruction of the KAM 
curves. In this case, we found a first order resonant island surrounded by a chaotic sea (see 
Fig. E](c)) and the reversal of the travel direction, for initial conditions falling outside the 
librational region, is now possible (we refer to this situation as bidirectional mixed chaos). 
As the ripple amplitude is increased further, the size of the region of librational motion 
becomes smaller, as shown in Fig. [2]^d). 

For sufficiently small ripple amplitudes, a particle executing librational motion, or ro- 
tational motion whose corresponding orbit in the PS lies near the separatrix, will collide 
with the walls almost perpendicularly, as Fig. |2]^a) shows. Under this circumstance, we 
can estimate the critical angle 9$ — @c for which the particle executes the largest librational 
motion. In this case, the reversal of travel direction occurs at a point that is arbitrarily close 
to the nearest hyperbolic fixed point x^ to the right of Xq if cos# c > 0. Since the average 
width of the channel is small, 6 C is virtually independent of the initial condition y (this is 
corroborated by our numerical results), and we can set the departure point at (x , 2/2(^0)) 
to estimate 8 C . Therefore, (xo,cos6 c ) is an initial condition in the PS. In order to estimate 
9 C , the condition of adiabatic invariance is used 

D(x n ) \sin8 n \ = D(x m ) |sin0 m | , (2) 

written for two arbitrary points (x n ,8 n ), (x m ,8 m ), with D(x n ) as the separation between 
the rippled walls at the point of collision x n . Setting x n = x , 6 n = 9 C , then x m ~ x^ and 
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Figure 2. (Color on line) Poincare plots for r = 1/2 and b = 0.1. (a) a = 0.001, (b) a = 0.02, (c) 
a = 0.04, and (d) a = 0.09. 

&m ~ § ■ Upon expanding through a Taylor series expression ([2]), but maintaining the terms 
of order a/b and (3%, where (3 C = n/2 — 9 C , we obtain 

(/3 C ) 2 w ^ [sin (2vr(xJ + r)) + sin (2ttx ) - sin (2-kx)) - sin (2tt(x + r))] . (3) 

Using the conservation of the energy and taking into account the condition of adiabatic 
invariance, which implies that the quantity C(x n ) = VoD(x n )\ sin^l remains constant, and 
with v representing the speed of the particle, we obtain 

l -( x f = E-V{x), (4) 

where E is the energy of the particle (with unit mass) and V(x) = \ ^S^) ' wn i cn can be 
interpreted as an effective potential that explains why the Poincare plots for a small ripple 
amplitude resemble the phase space of a simple pendulum. The hyperbolic fixed points can 
be computed by obtaining the maxima of V(x). In our specific model (r = 1/2), x h j = 3/4 
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if xq = and cos# c > 0. In this case, expression fl3]) is reduced to 



b J 

The condition of adiabatic invariance and the critical angles have been used previously 
to arrive at estimates of the electronic transport properties in rippled channels when the 



Poincare plots of the associate dynamical system exhibit a pendulum-like phase space 23 



24]. 



III. THERMAL PROPERTIES OF THE CHANNEL AT SMALL RIPPLE AMPLI- 
TUDE 

This section considers a finite version of the channel discussed in the previous section. 
This channel consists of N replicas (a replica is a fundamental cell formed with the rippled 
walls and the vertical dashed lines, as shown in Fig. [1]) of length 1 = 1. One end of the 
channel is at xl = 0, so the other one is at xr = N. To induce heat transport, the channel 
is placed between two heat baths at temperatures Tl and Tr (Tl > Tr, see Fig. [[]), modeled 
by stochastic kernels of Gaussian type by 

PM - (-|) , P W = -J= exp (-|) , (6) 

where P(v x )P(v y ) is the probability distribution of the velocities for the particles emerging 
from the baths. The values of N are such that N ^> 6; hence, the number of particles that 
can cross the channel without colliding with the rippled walls is much lower than that of 
those which will have collisions. 

Due to the fact that the particles do not interact among themselves, the motion of a 
single particle is considered over a long time period as it collides with the rippled walls and 
thermal baths. Collisions with the rippled boundaries are specular, and the velocity of the 
particle just after a collision with a thermal bath is determined by the distribution given by 



Eq. (jfJJ). Following the ideas from Ref. J2], we divide the configuration space into slices {Cj} 
(eventually a slice will be taken as a fundamental cell). The time that the particle spends 
in the slice in the j-th visit is denoted by tj, and the total number of times that the particle 
crosses C{ is M. The temperature of the slice Cj is defined by 



where Ej(Ci) is the kinetic energy of the particle at the j-th visit of the slice Cj. The heat 
flux of one particle can be defined by 

(8) 



1 Nc 



k=l 



where (AE) k is the energy change of the particle at the k-th collision with a thermal bath, 
and t c is the total time that the particle takes to collision N c times with a thermal bath. 
The correct way to obtain the thermodynamic limit is to keep the number of particles per 
cell fixed as the size N of the channel increases. For instance, if a single particle per cell is 
considered, then the heat flux in a channel of N length is jAr(t c ) = Nji(t c ). The thermal 
conductivity k is determined by 

N 2 h(N) 



K = 



N 2 ji(N), 



(9) 



T L -T R 

where we have defined the thermal conductivity through the Fourier law of heat conduction. 
The next step is to obtain an expression of the temperature profile and heat flux in the 
stationary regime for the case of a small ripple amplitude. In the interests of simplicity, but 
without losing generality, the slices {Ci} are assumed to be equal to the fundamental cells. 
For sufficiently small ripple amplitudes, the librational orbits, or rotational orbits near the 
separatrix, represent a particle bouncing off the walls almost perpendicularly. Thus, the 
time in the j-th visit that the particle spends within the slide C\ (which is in direct contact 
with the hot bath) for these rotational orbits is given, approximately, by 

dx 



t 



(r) 



E 

Ci 

i 



JJ) JJ) 
x k+l x k 
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Ji) 
x k 


Jo 



I 



D(x)dx 



.0") 



D 2 (x) -D 2 (0) sm 2 (0 { j>) 



(10) 



while for the librational orbits, we have 



Ci 

= 2 / 

Jo 



Si) 



X k+1 X k 



(J) 



d3) 



X 



dx 

l ^Ej-V^x) 
D(x)dx 



(11) 



otfyD 2 (x)-D 2 (0)fan 2 (e®) 

x k\ x i^) denotes the velocity in the x-direction just after the k-th collision with any wall (the 
re-coordinate of the particle at the k-th collision with any wall) during the j-th visit to cell C\. 
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function of xjj? is determined using the effective potential given by Eq. (j3J), while 
the trajectory of the particle becomes known once the random initial conditions Vq,9q 
are determined just after a collision with a thermal bath. When the particle collides with 
a thermal bath again, we continue computing its trajectory once the new initial conditions 
Vq , 6q are known. Here 6q\ Vq \ Ej, Vj(x) have the same meaning as #0, vq, E, V(x), 
introduced above, respectively. D(x) is the distance between the rippled walls at point x, 
which in our specific model is given by 



Expressions ( flQj) and ( fTTl) are obtained using Eq. (jl]), where we set C® = D(0)v^'\sme^\. 
The factor 2 in Eq. (flTl) derives from the fact that the time that the particle takes to 
execute librational motion before being reinjected into the same bath is twice the time it 
spends going from the thermal bath to the returning point x c , where the particle changes 
its direction of travel and the corresponding x-velocity is approximately zero. Then, x c is 
determined by the nearest positive root of E — V(x) = to x = 0. In the case of the 
slice C^, we find similar expressions, but x c is now the nearest root to the left of xq = N 
(note that, due to the system periodicity, this is equivalent to take the absolute value of 
the nearest negative root to Xo = 0) that we denoted as x c n- Particles emerging from the 
thermal baths cannot have access to librational orbits in Cj cells with % = 2, 3, . . . iV — 1; 
thus the motion in these cells is rotational. 

Rotational orbits that do not lie near the separatrix are almost flat as shown in Fig 
2(a); so ip ~ (vq cosO^)^ 1 . A rotational orbit near the separatrix means that its corre- 
sponding initial condition satisfies cos# c ~ (a /2b) 1 < cosOq^ . For rotational orbits that 
do not lie near the separatrix, tj can also be approximated by Eq. ( ITUi) . since upon de- 
veloping it in power series of a/b, we obtain tj = (v^ cosO^)^ 1 + O (^/(tfcos^O^)) + 



O [a 2 / (b 2 cos 2 9q ) ) . Due to the periodicity of rotational motion, the time that the particle 
takes to cross any cell while executing librational motion, is approximated by Eq. ( 110)) . 



Equation fllOp can be developed in power series of (a/t&cos 2 ^) + a 2 /(6 2 cos 2 ^ ) )J, but 
we decide not to present the expansion because as we approached an initial condition that 
falls near the separatrix, it became necessary to include more terms of the expansion; while 
for an initial condition that falls exactly on the separatrix all terms of the expansion must 
be included. When we consider librational motion, another expansion is necessary; so we 



D(x) = 2(6-asin(27rx)). 
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found it more convenient to keep Eqs. (fTUj) and (llip in their integral form, 
ow 

and tj, which is given by 



(r) 

Now we can establish the behavior of tj as a function of both the cell number and of tj 



7(1 

if, if 9 >(), i I. A" 



where C is the critical angle defined in Eq. (j5]). Expression f fl2|) is derived using the 
arguments that follow. The particles that come from the thermal baths at an angle 9 > 6 C 
execute librational motion in the cells i = 1 or i = N, depending on the thermal bath from 
which they come. Those particles are then reinjected into the same heat reservoir, then tj 
is only defined in the cells i — 1, N. 9 C is the same for particles coming from any heat bath 
because the spatial separation of the baths is equal to an integer number. When 9 < 9 C , 
particle motion is rotational and therefore periodic; thus the time required for a particle to 
cross cell Ci is independent of the i index. 

To compute the temperature profile, we divide the particles in two types: those that 
come from the bath at temperature Tl and those that come from the bath at temperature 
Tr at the moment they cross the slice Ci. In this way, definition (JTj) reads 



T(d) = ' ,' R , (13) 

<*>? + <*>i? 



(tE)f + (tE) ( S 

± \yi) - 

where (• • • )® denotes the time average at cell Ci using probability distribution ([6]) with 
T = T L ; (■■■ )® has the same notation, but with T = Tr. By expressing the probability 
distribution (J5]) in polar coordinates (v ,9 ), an d substituting Eq. ( 1T2"j) into Eq. ( TT3"|) . we 
have 



T(C t )~^/nTn~ for i = 2, 3, . . . , N - 1, (14) 

T{Cl) * VTLTR VT- R h(e c ) + (VT L + VTdU9 c y (15) 

T(C N ) * VT L T R ^ Lm) + {VT - +VT - )hidcY (16) 



Here, 



p c rx c (M,(x cN (M) D(x)sm(p )dxdp 



^/D 2 {x)-D 2 {0) cos 2 (/3 ) 
D{x) cos(9 )dxd9 



o y/D*(x)-D 2 (0)sm 2 (6 ) 
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with [3q = | — 6q. The numerical results support our theoretical results for the temperature 
profile predicted by Eqs. f|T4|) - f|T6|) . as shown in Fig. [3]^a), where it is clear that the tem- 
perature of the cells that are not in direct contact with the thermal baths is the geometric 
average of the temperatures of those baths. Deviations from this flat profile occur at the 
cells that are in direct contact with the heat reservoirs due to the librational motion of the 
particles in those cells. The theoretical predictions for these deviations agree with the nu- 
merical simulation. In the case of a flat channel (a = 0) 9 C — f and the temperature profile 
is completely flat with value a/T^Tr, this result is exact and can be corroborated directly 
using definition ( [131) . 

Using a similar procedure that made it possible to obtain analytical results for the tem- 
perature profile, we obtain the following expression for heat flux 



where the resulting factor N in the denominator comes from the periodicity of the rotational 
motion, and from the fact that a particle executing this kind of motion must travel through 
N cells of unit length before colliding with the other thermal bath. At the limit N ^> 1, 
heat flux becomes independent of system size and thermal conductivity scales as k ~ N. 
We have corroborated numerically that the latter scaling behavior is valid as long as the 
system exhibits regular dynamics. In the latter case, ballistic diffusion is also present. 

Our theoretical prediction ( TT7l) is confirmed by the numerical results, as the inset in Fig. 
Et^b) shows. We should point out an unexpected behavior of increasing heat flux when a is 
increased to a maximum value at a »s 0.015, when chaotic behavior appears. This behavior 
can be explained as follows: the particles executing librational motion do not contribute to 
heat flux because they never reach the other side of the channel. This fact is associated 
with the term sin0 c in the numerator of Eq. (fT7|) . which decreases as ripple amplitude is 
increased. However, the average time that the particle takes to travel from one thermal 
bath to the other becomes smaller as ripple amplitude increases, because only particles with 
larger velocity components in the x-direction can reach the other side of the channel. The 
reduction of the average time is associated with the term J 2 in the denominator of Eq. (ITT)) ; 
thus, we have two competing mechanisms and the strongest one is that associated with the 



3n(N) 



m L +m R +N((&)) L+ (&)) R ) 




(17) 



2§ (V^7i(0 c ) + y/TEl 3 (9 c ) + N (VTE + V7r) Wc)) ' 
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term I2 when N ^> 1. 
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Figure 3. (Color on line) (a) Temperature profile for three different amplitudes; a = 0.001, a = 
0.003 and a = 0.01. Continuous lines represent the theoretical prediction ()14|) -f)16[) . while the 
different styles of dots represent the numerical data. The inset presents a zoom of the numerical 
temperature profile in the central part of the channel to show the accuracy of the theoretical 
prediction (|14p . (b) Heat flux of one particle as a function of ripple amplitude. The inset uses 
dots to show the numerical data, while the continuous line represents the theoretical prediction for 
one particle heat flux for the amplitude interval [0,0.01]. In both figures, Tl = 10, Tr = 1, and 
channel length is N = 10. 



The behavior of the temperature profile predicted by Eqs. (j!4p - (jl6p is similar to the 
corresponding one for a homogeneous harmonic chain between two Langevin heat baths 



25] ; while the scaling behavior of the thermal conductivity of the homogeneous chain and 
our model is identical. The main difference between the two systems comes from the fact that 
the bulk temperature of the harmonic chain is the arithmetical average of the temperatures 
of the two heat reservoirs. However, in both systems, deviations from a flat profile occur in 
those parts of the systems that are in direct contact with the reservoirs. Another system 
that presents the same size-scaling behavior in heat flux and thermal conductivity is the 
harmonic chain, with specific long-range correlation of the isotopic disorder, attached to 



Langevin heat baths |2f 
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IV. LARGER RIPPLE AMPLITUDE 



We consider four channels, with different values of a that produce different degrees of 
chaos in the system, in order to study their thermal transport properties. In case I, a — 0.02 
and unidirectional mixed chaos is present (see Fig. EJ^b)). Case II corresponds to a = 0.04, 
where the system exhibits bidirectional mixed chaos (Fig. EJ^c)). Cases III and IV correspond 
to a = 0.07 (not shown) and a = 0.09 (Fig. EJ^d)), respectively; bidirectional mixed chaos 
is presented as well, but with a smaller region of librational motion due to stronger chaotic 
behavior. 

We are interested in studying two different, but related quantities: first, the scaling 
behavior of thermal conductivity in relation to the system size; and, second, the mean square 
displacement (Ax 2 ) = {(x(t) — a;(0)) 2 ), where (...) denotes an average over different initial 
conditions and x(t) is the x-coordinate of the particle position at time t. For our numerical 
simulations of heat conduction, we choose Tl = 10 and Tr = 1. It is important to note 
that the scaling of ji(N) with system size does not depend on the properties of the thermal 
baths, but exclusively on the geometry of the channel [6j. Thus, the thermal conductivity 
defined by Eq. gives the same scaling behavior regardless of the temperature difference 
between the reservoirs. At the channel lengths explored in the numerical simulation, the 
thermal conductivity can be approximated by k = AN 13 (A and constants) for sufficiently 
long channels (see Fig. H^a)). In the case of unidirectional mixed chaos, the exponent is 



close to 1; ballistic behavior as reported in Ref. 22] is not present due to the presence of 
particles whose initial conditions fall in the chaotic separatrix and therefore their motion is 
not unidirectional. As soon as bidirectional mixed chaos appears, an abrupt decay occurs 
in exponent from « 1 to ~ 0, because there are no more particles with unidirectional 
motion. As ripple amplitude is increased, the exponent approaches zero. The diffusion 
of particles in the corresponding channel of infinite length is also studied, and within the 
time interval explored in Fig. Hl^b), the growth of the mean square displacement can be 
approximated by Dt a , for a sufficiently long time period. In reference to the scaling behavior 
of the thermal conductivity and mean square displacement, we say "can be approximated", 
and not "scale as", because exponent exhibits a small change in relation to system size, 
and exponent a presents a small change in relation to time t. This means that these two 
exponents remain practically constant for relatively long channel interval lengths and long 
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time periods, respectively. The origin of the size- dependent and time-dependent exponents 
is discussed in the following paragraph. The numerical simulations of particle diffusion 
concord with the corresponding simulation of heat flux, in that exponents (3 and a (see Figs. 
Ufa) and|U(b)) satisfy the relation (3 = 2 — 2/ a, derived in Ref. [5(. It is interesting to note 
that for bidirectional mixed chaos, exponent a is close to 1 and exponent is close to 0. 




Figure 4. (a) (Color on line) Heat flux of one particle as a function of system size N in log-log 
scale, (b) Mean square displacement (Ax 2 (t)^ against time t. 10 5 particles are used to obtain 
the diffusive properties of the system. Particles are initially at x = 0, initial velocities obey the 
Maxwell-Boltzmann distribution at temperature T = 1. In both figures, case I is (x), case II (□), 
case III (O) and case IV (A). 

In any model with an infinite horizon, there are particles that can travel arbitrary dis- 
tances without suffering any collision with the boundaries at time t. These particles make a 
contribution to the velocity auto-correlation function that scales as 1/t, and is the leading 
term if the corresponding contribution of the rest of the particles (z.e., those that do collide 
at time t) decays faster than 1/t 2^]. The algebraic decay 1/t of the velocity auto-correlation 
function has been confirmed numerically for the Lorentz gas with an infinite horizon 28- 
301 ]. and for a polygonal channel with an infinite horizon [lfj. Therefore, in our model, the 
mean square displacement must grow as ct log (t) + dt (the diffusion constant is related to 
the velocity auto-correlation function by the Einstein-Green-Kubo formula) if the degree 
of chaos is strong enough; such a correction will appear in Fig. ID^b) as a small change in 
the a exponent and, therefore, a curvature in the mean square deviation will be visible for 
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longer times. To observe the logarithmic correction, we plot ((Ax 2 )/ 1) against log(t), as 
shown in Fig. [51 and found that for cases II, III and IV, (Ax 2 ) = ct log (t) + dt with < c 
and d ^ c (in general, we find this behavior for bidirectional mixed chaos). For case I, the 
scaling law ct log (t) + dt is incorrect; in fact, we do not know the correct asymptotic scaling 
law for unidirectional mixed chaos. However, (Ax 2 ) = Dt a (with a and D constants) is a 
good approximation at the time period explored in the numerical simulations. The ratio 
c/d decreases as ripple amplitude is increased, this is simply a consequence of the reduction 
of the number of particles that can travel arbitrarily far without colliding with the walls, as 
ripple amplitude becomes larger. In this sense, we can say that diffusion approaches normal 
behavior as ripple amplitude increases; hence, the same statement is also valid for heat flux. 
The marginally anomalous diffusion exhibited by our system (in the case of bidirectional 
chaos) provides evidence that the contribution to the velocity auto-correlation function from 
particles having collisions at time t decays faster than 1/t. Then, if the particles that can 
travel arbitrarily far without interacting with the walls were not present, diffusion should 
be normal. This argument indicates that normal diffusion can take place in systems with 
mixed phase space when there are not KAM curves that preclude connection of all chaotic 
regions, and when non-chaotic trajectories do not contribute to heat transport. 

It must be stressed that for cases III and IV, all initial conditions of the particles coming 
from the thermal baths fall on the chaotic sea. However, if we place the thermal baths such 
that the initial conditions of these particles can have access to periodic and quasi-periodic 
orbits, exponents /3 and a should be insensitive to this change because particles executing 
periodic or quasi-periodic motion do not contribute to heat flux for they are reinjected 
into the same bath from which they originally came. We have tested the latter statement 
numerically by placing the thermal baths at xl = 0.25 and xr = 0.25 + N. Here, particles 
have access to regular motion for all different cases. We find that exponents /3 and a are 
practically the same than the corresponding ones given by Fig. HI We report the following 
values for exponent (3: case I, j3 — 0.995; case II, (3 = 0.089; case III, (3 = 0.044, and case 
IV, = 0.023. 

The behavior of the temperature profile for four different cases is shown in Fig. [6j It 
is clear that in case I the temperature profiles in the central part of the channel are linear 
and maintain a similar shape for different system sizes. This behavior is quite similar to 
that presented in references jf], |2l| , when the first system have no disorder and the second 
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Figure 5. (Color on line) f(t) = ((Ax 2 ) /t) as a function of log(i). Figures from top to bottom 
correspond to cases II, III and IV, respectively. Solid lines represent the best fit of the numerical 
data. The fit gives the following values: case II, c = 0.18 ± 1.9 x 1(T 2 , d = 10.00 ± 0.2 ; case III, 
c = 9.24 x 10" 3 ± 1.8 x 10~ 3 , d = 0.67 ± 1.4 x 10~ 2 ; case IV, c = 6.50 x 10~ 3 ± 3.47 x 10~ 4 , 
d = 8.45 x 10~ 2 ± 1.3 x 10~ 4 . In cases II and III, the statistical errors are less than the size of the 
points; in case IV, statistical errors are represented by the error bars. 

one exhibits regular dynamics. There are jump discontinuities in the temperature profile 
for cases I and II at the ends of the channel. This phenomenon is the result of boundary 
heat resistance that usually appears when there is a heat flux across the interface of two 
adjacent materials (see. Ref. 3l| and references therein). These jumps are finite size 
effects because as the system size increases boundary heat resistance decreases, as do the 
temperature jumps. This is seen clearly in Fig. [6(b). When /3 ~ there are almost no 
jumps (see Figs. Etc ) andE(d)) and the temperature profile can be approximated by the 



formula (see Ref. 



2l| for details): 



T{x) 



T L T£(1 



(1 -x)t + T l 2 ^ 



(18) 



3 

which is valid for large system sizes. Here x = j?, 7 = (2 — a) a 2 and a is the diffusion 
exponent (keeping in mind that the diffusion exponent is not a constant, though it can be 
treated as such for relative long periods of time). For case I, the temperature profile predicted 
by (|T8|) is clearly different from our numerical simulation. The numerical simulation and 
temperature profile predicted by this formula have a similar shape but are displaced (see 
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Eq. (18) 
K=64 
N=128 
N=512 




Figure 6. (Color on line) Internal local temperature as a function of the rescaled cell number i/N 
with Tl = 10 and Tr = 1. (a) Corresponds to case I; (b) to case II; (c) to case IV; and, (d) also 
to case IV, but with the thermal baths placed at xl = 0.25 and xr = 0.25 + N. The solid lines 
in panels (c) and (d) correspond to the best fits for the numerical temperature profile at N = 192 
with Eq. p8p. the fit results are the following: panel (c), a = 1.004; panel (d), a = 1.005 . In 
panel (a), the solid line is the plot of Eq. (|18p with a = 1.99 and the dashed line represents the 
best fit for the temperature profile at N = 512 with the heuristic formula (|19p . the result of the fit 
is a = 1.99. In panel (b), we set a = 1.00 and it is clear that the temperature profile approaches 
the asymptotic temperature profile (|18p as the length of the channel is increased. 



Fig. E(a)). However, we can see that the heuristic formula 



T(x) 



T L T£(l-x)i + T R Tlx^ 



THi-x^ + Tlx-y 



(19) 



shows good agreement with the numerical results for case I (see dashed line in Fig. [6^a)). 
Even in the case of the flat channel (a = 2), there is a discrepancy between formula f|T8|) 
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and (I19p . since the former predicts a flat profile given by T(x) = fggf-, while formula ([19]) 
also predicts a flat temperature profile with the value T(x) = a/Tr7Y, but the latter flat 
temperature profile is an analytical exact result that can be obtained directly from definition 
(J7j). When the formula ( |T8l) was derived, a definition of temperature in terms of the mean 
density of particles at the cell was used. Nevertheless, definition ([7]) is the time average 
of kinetic energy at cell Cf, therefore, these two definitions lead to different temperature 
profiles that are only equal when diffusion is normal. This explains the clear discrepancy 
between our numerical simulations and Eq. (1T81) for case I. In any case, both definitions 
of temperature lead to temperature profiles that are closely related to diffusion exponent 
a for sufficiently large system sizes. If a system obeys Fourier's law (a = 1), then the 
formula derived in Ref. [2| is recovered. In addition, when the temperature difference of the 
reservoirs is small, the typical linear temperature profile will be obtained. 



V. CONCLUSIONS 

We analyzed the thermal transport properties of a sinusoidal channel placed between 
two Gaussian type thermal baths. The dynamics of a particle moving in the corresponding 
infinite length channel exhibits Poincare plots with a KAM structure typical of generic 
Hamiltonian systems. When the ripple amplitude is small, the dynamics of the system is 
regular. For this case, estimates of the heat flux and temperature profile were obtained using 
an effective potential. Specifically, the temperature profile of the central part of the channel 
is well approximated by a flat profile given by the geometric average of the temperatures of 
the two reservoirs. Consequently, the same result is to be expected for other billiard systems 
with similar dynamics. In the regime of regular dynamics, the thermal conductivity scales 
with the system size as k ~ iV for sufficiently large system size and the diffusion of particles 
is ballistic. Unidirectional mixed chaos appears as the ripple amplitude is increased; and 
the mean square displacement can be approximated by Dt a (being a < 2 a constant) for 
a relative long time period. When bidirectional mixed chaos appears, the mean square 
displacement grows asymptotically as t log (i), then diffusion exponent is time dependent, 
but it remains practically constant for relative long time periods. Temperature profiles were 
also analyzed for different degrees of chaos in the system, and it was found that the diffusion 
exponent is closely related to the temperature profile of the system. 
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